#include "Lattice.hpp"
#include <iostream>
#include <map>
#include <fstream>

Lattice::Lattice(uint lengthx,uint lengthy,real &omega):
	lengthx_(lengthx),
	lengthy_(lengthy),
	omega_(omega)
{
	std::cout<<"lengthx: "<<lengthx_<<"\n";
	std::cout<<"lengthy: "<<lengthy_<<"\n";
	dest_ = PDF_Field(lengthx_,lengthy_);
	src_ = PDF_Field(lengthx_,lengthy_);
	vs_ = V_Field(lengthx_,lengthy_);
	rhos_ = D_Field(lengthx_,lengthy_);
	nodesType_ = Flags(lengthx_,lengthy_);

	for(uint i=0;i<lengthx_;++i)
	{
		for(uint j=0;j<lengthy_;++j)
			nodesType_(i,j)=FLUID;
	}

}

Lattice::Lattice(real &omega, std::string &pgm_file):
	omega_(omega)
{
	std::ifstream matrix_reader(pgm_file);
	if (matrix_reader.is_open())
	{
		std::string name;
		int max_value;
		matrix_reader >> name >> lengthx_ >> lengthy_ >> max_value;
		std::cout<<"lengthx: "<<lengthx_<<"\n";
		std::cout<<"lengthy: "<<lengthy_<<"\n";
		dest_ = PDF_Field(lengthx_,lengthy_);
		src_ = PDF_Field(lengthx_,lengthy_);
		vs_ = V_Field(lengthx_,lengthy_);
		rhos_ = D_Field(lengthx_,lengthy_);
		nodesType_ = Flags(lengthx_,lengthy_);

		// Set the Obstacle as per the Nodes
		for (int j = lengthy_-1; j >=0; --j) 
		{
			for (int i = 0; i < lengthx_; ++i) 
			{
				int k;
				matrix_reader >> k;
				if (k/max_value==1) nodesType_(i,j) = FLUID;
				else nodesType_(i,j) = OBSTACLE;
			}
		}
	}
	matrix_reader.close();
}





void Lattice::initialize(real density, real velX, real velY)
{
	inversePDF.insert(std::pair<uint,uint>(NORTH,SOUTH));
	inversePDF.insert(std::pair<uint,uint>(SOUTH,NORTH));
	inversePDF.insert(std::pair<uint,uint>(EAST,WEST));
	inversePDF.insert(std::pair<uint,uint>(WEST,EAST));
	inversePDF.insert(std::pair<uint,uint>(NWEST,SEAST));
	inversePDF.insert(std::pair<uint,uint>(SEAST,NWEST));
	inversePDF.insert(std::pair<uint,uint>(SWEST,NEAST));
	inversePDF.insert(std::pair<uint,uint>(NEAST,SWEST));
	inversePDF.insert(std::pair<uint,uint>(CENTER,CENTER));

	weights.insert(std::pair<uint,real>(CENTER,4.0/9.0));
	weights.insert(std::pair<uint,real>(EAST,1.0/9.0));
	weights.insert(std::pair<uint,real>(NEAST,1.0/36.0));
	weights.insert(std::pair<uint,real>(NORTH,1.0/9.0));
	weights.insert(std::pair<uint,real>(NWEST,1.0/36.0));
	weights.insert(std::pair<uint,real>(WEST,1.0/9.0));
	weights.insert(std::pair<uint,real>(SWEST,1.0/36.0));
	weights.insert(std::pair<uint,real>(SOUTH,1.0/9.0));
	weights.insert(std::pair<uint,real>(SEAST,1.0/36.0));


	for(uint i=0;i<lengthx_;++i)
	{
		nodesType_(i,0)=OBSTACLE;
		nodesType_(i,lengthy_-1)=OBSTACLE;
	}
	for(uint j=0;j<lengthy_;++j)
	{
		nodesType_(0,j)=OBSTACLE;
		nodesType_(lengthx_-1,j)=OBSTACLE;
	}

	real tmpFe=0;
	for(uint i=0;i<lengthx_;++i)
	{
		for(uint j=0;j<lengthy_;++j)
		{
			if(isFluid(i,j))
			{
				vs_(i,j,0)=velX;
				vs_(i,j,1)=velY;
				rhos_(i,j)=density;
				nodesType_(i,j)=FLUID;
				for(uint k=0;k<DIMEN;++k)
				{
					tmpFe=fequilibrium(i,j,k);
					dest_(i,j,k)=tmpFe;
					src_(i,j,k)=tmpFe;
				}
			}
		}
	}
}

void Lattice::streaming(real vwx,real vwy)
{
	uint x=0;
	uint y=0;
	for(uint i=1;i<lengthx_-1;++i)
	{
		for(uint j=1;j<lengthy_-1;++j)
		{
			if(isFluid(i,j))
			{
				for(uint k=0;k<DIMEN;++k)
				{
					x = i + eX[k];
					y = j + eY[k];
					if(isFluid(x,y))
						dest_(x,y,k)=src_(i,j,k);
					else
					{
						if(y==lengthy_-1)
							dest_(i,j,inversePDF[k])=src_(i,j,k)-6*weights[k]*(vwx*eX[k]+vwy*eY[k]);
						else
							dest_(i,j,inversePDF[k])=src_(i,j,k);
					}                

				}

			}
		}   
	}

}

void Lattice::collision()
{
	real tmpFe=0;
	for(uint i=1;i<lengthx_-1;++i)
	{
		for(uint j=1;j<lengthy_-1;++j)
		{
			if(isFluid(i,j))
			{
				updateDensity(i,j);
				updateVelocity(i,j);
				for(uint k=0;k<DIMEN;++k)
				{
					tmpFe=fequilibrium(i,j,k);
					dest_(i,j,k)=dest_(i,j,k)-omega_*(dest_(i,j,k)-tmpFe);
				}
			}   
		}
	}
}

void Lattice::streamingCollision(real vwx,real vwy)
{
	real tmpFe=0;
	uint x=0;
	uint y=0;
	for(uint j=1;j<lengthy_-1;++j)
	{
		for(uint i=1;i<lengthx_-1;++i)
		{
			if(isFluid(i,j))
			{
				for(uint k=0;k<DIMEN;++k)
				{
					x = i + eX[k];
					y = j + eY[k];
					if(isFluid(x,y))
						dest_(i,j,k)=src_(x,y,k);
					else
					{
						if(y==lengthy_-1)
							dest_(i,j,k)=src_(i,j,inversePDF[k])-6*weights[k]*(vwx*eX[k]+vwy*eY[k]);
						else
							dest_(i,j,k)=src_(i,j,inversePDF[k]);
					}
				}

				updateDensity(i,j);
				updateVelocity(i,j);
				for(uint k=0;k<DIMEN;++k)
				{
					tmpFe=fequilibrium(i,j,k);
					dest_(i,j,k)=dest_(i,j,k)-omega_*(dest_(i,j,k)-tmpFe);
				}
			}
		}
	}
}

void Lattice::copyGrid()
{
	for(uint i=1;i<lengthx_-1;++i)
	{
		for(uint j=1;j<lengthy_-1;++j)
		{
			if(isFluid(i,j))
			{       
				for(uint k=0;k<DIMEN;++k)               
					src_(i,j,k)=dest_(i,j,k);
			}
		}
	}
}

void Lattice::comparePDFsGrids()
{
	real sum=0;
	uint counter=0;
	for(uint i=1;i<lengthx_-1;++i)
	{
		for(uint j=1;j<lengthy_-1;++j)
		{
			if(isFluid(i,j))
			{       
				for(uint k=0;k<DIMEN;++k)               
					sum+=src_(i,j,k)-dest_(i,j,k);
				counter+=DIMEN;
			}
		}
	}
	std::cout<<"Difference: "<<sum/counter<<"\n";
}

void Lattice::printF(const std::string &vtkFile)
{
	VTKFileWriter vtk(lengthx_,lengthy_,vtkFile);
	vtk.WriteVector(vs_,"Velocity");
	vtk.WriteScalar(rhos_,"Density");
	vtk.WriteFlag(nodesType_,"FLAGS");
}

void Lattice::checkMaps()
{
	for(uint k=0;k<DIMEN;++k)
	{
		std::cout<<"direction "<<k<<" inverse " <<inversePDF[k]<<" weight "<<weights[k]
		<< " unitiyVec ("<<eX[k]<<","<<eY[k]<<")\n\n";
	}
}

